module fonctions
implicit none
!------------------------------------------------------------------------
double precision, parameter :: x = 20.d0 !Mesh refinement: x smaller => mesh smaller without syst.size modification
double precision, parameter :: y = 0.7d0 !System size: y bigger => system bigger without mesh modification
!Adjust spatial mesh to prevent instabilities (fine enough), reasonable CPU time (not too fine), and boundaries
!Euler stability criterion: (DX)^2/(D*(Dt)) > 2, D=10^{-9} [m^2/s], DX=realspace, Dt=realtime
!------------------------------------------------------------------------
double precision, parameter :: U = -38.64d0*y !Imposed potential: Nl*y => U*y
double precision, parameter :: deltat = 0.1d0 !time discretisation.
integer, parameter :: Nsavedata = 200 !Saving data every Nsavedata iterations
integer, parameter :: nsave = 2 !maximum number saves. Then last one overwritten. Alonissos: =2. Amorgos: =15.
!------------------------------------------------------------------------
double precision, parameter :: realspace = .01d0*x !Conversion factor for real space [centimeter]: /x
double precision, parameter :: realtime = (10.d0/60.d0)*x**2 !Conversion factor for real time [minutes]: /x^2
integer, parameter :: Nl = nint((10000.d0/x)*y) !#sites to left; mesh size: Nl*realspace=constant | /x
integer, parameter :: Nr = nint((10000.d0/x)*y) !#sites to right; mesh size: keep Nr*realspace=constant | /x
!------------------------------------------------------------------------
double precision, parameter :: a0 = 10.d0 !Initial concentration of A
double precision, parameter :: b0 = 0.1d0 !Initial concentration of B
double precision, dimension(1:4), parameter :: D = (/ 0.1d0, 1.d0, 1.d0, 1.d0 /) !Diffusion constants for A+,A-,B+,B-
!------------------------------------------------------------------------
logical, parameter :: infiniterate = .true. !rate of reaction: infinite or finite
double precision, parameter :: rate = 100.d0 !rate of reaction if finite
!------------------------------------------------------------------------
character*(*), parameter :: filename0 = 'liese-state.txt'
character*(*), parameter :: filename1 = 'liese-C.txt' !concentration of C(x,t): t,{xj,C(xj)}
character*(*), parameter :: filename2 = 'liese-R.txt' !reaction rate R(x,t): t,{xj,R(xj)}
character*(*), parameter :: filename3 = 'liese-I.txt' !current I(t): t,I(t)
character*(*), parameter :: filename4 = 'liese-E.txt' !electric field E(t): t,{xj,E(xj)}
character*(*), parameter :: filename4b = 'liese-Et.txt' !electric field E(t) at the front: t,E(xfront)
character*(*), parameter :: filename5 = 'liese-ap.txt' !concentration of A+: t,{xj,A+(xj)}
character*(*), parameter :: filename5b = 'liese-apt.txt' !concentration of A+ at the front: t,A+(xfront)
character*(*), parameter :: filename6 = 'liese-am.txt' !concentration of A-: t,{xj,A-(xj)}
character*(*), parameter :: filename7 = 'liese-bp.txt' !concentration of B+: t,{xj,B+(xj)}
character*(*), parameter :: filename8 = 'liese-bm.txt' !concentration of B-: t,{xj,B-(xj)}
character*(*), parameter :: filename8c = 'liese-bmt.txt' !concentration of B- at the front: t,B-(xfront)
character*(*), parameter :: filename8b = 'liese-AB.txt' !concentration of all species-: t,{xj,all(xj)}
character*(*), parameter :: filename9 = 'liese-q.txt' !charge: t,{xj,(A+ - A- + B+ - B-)(xj)}
character*(*), parameter :: filename10 = 'liese-Qtot.txt' !total charge Q(t): t,Q(t)
character*(*), parameter :: filename11 = 'liese-Cfront.txt' !position of the C front: t,Cfront(t)
double precision, parameter :: Cfront = 0.0001 !If C < Cfront then this is the position of the C front
integer, parameter :: Ns = Nl + Nr !total number of sites. Do not adjust.
real, dimension(1:nsave,1:Ns) :: tmpC,tmpE,tmpR,tmpAP,tmpAM,tmpBP,tmpBM,tmpQ,tmpAB !temporary for data saving.
integer :: pouet
contains

!------------------------------------------------------------------------
subroutine DerivativesRHS(AB,E,Dxx,Ex,XI)
double precision, dimension(1:4,1:Ns), intent(in) :: AB
double precision, dimension(1:Ns), intent(inout) :: E
double precision, dimension(1:4,1:Ns), intent(inout) :: Ex,Dxx
double precision, intent(out) :: XI
double precision, dimension(1:Ns) :: sum1,sum2 !Sums in expr. of E
double precision :: tmp
integer :: i
!E field
do i=2,Ns-1
	sum1(i) = D(1)*(AB(1,i+1)-AB(1,i)) - D(2)*(AB(2,i+1)-AB(2,i)) &
					+ D(3)*(AB(3,i+1)-AB(3,i)) - D(4)*(AB(4,i+1)-AB(4,i))
	sum2(i) = (D(1)*(AB(1,i+1)+AB(1,i)) + D(2)*(AB(2,i+1)+AB(2,i)) &
					+ D(3)*(AB(3,i+1)+AB(3,i)) + D(4)*(AB(4,i+1)+AB(4,i)))/2.d0
end do
sum1(1)=0.d0
sum2(1)=a0*(D(1)+D(2))
sum1(Ns)=0.d0
sum2(Ns)=b0*(D(3)+D(4))
tmp=0.d0
do i=1,Ns
	tmp = tmp + 1.d0/sum2(i)
end do
XI=U/tmp
E(:) = (XI+sum1(:))/sum2(:)
!Right-hand side: $d^2 c_i / d x^2$ (Dxx) and $-z_i d / dx c_i E$ (Ex)
do i=2,Ns-1
	Ex(1:4:2,i) = - (E(i)*(AB(1:4:2,i+1)+AB(1:4:2,i)) - E(i-1)*(AB(1:4:2,i)+AB(1:4:2,i-1)))/2.d0
	Ex(2:4:2,i) = + (E(i)*(AB(2:4:2,i+1)+AB(2:4:2,i)) - E(i-1)*(AB(2:4:2,i)+AB(2:4:2,i-1)))/2.d0
	Dxx(:,i) = AB(:,i+1) - 2.d0*AB(:,i) + AB(:,i-1)
end do
return
end subroutine DerivativesRHS
!------------------------------------------------------------------------
subroutine RK(AB,E,Dxx,Ex)
double precision, dimension(1:4,1:Ns), intent(inout) :: AB
double precision, dimension(1:Ns), intent(in) :: E
double precision, dimension(1:4,1:Ns), intent(in) :: Ex,Dxx
double precision, dimension(1:4,1:Ns) :: tmpAB,p1,p2,p3,p4,tmpEx,tmpDxx
double precision, dimension(1:Ns) :: tmpE
double precision :: tmp
integer :: i
	do i=1,4
		p1(i,:) = D(i)*(Dxx(i,:)+Ex(i,:))
	end do
	tmpAB = AB+p1*deltat/2.d0
	tmpAB(1,:) = tmpAB(2,:) + tmpAB(4,:) - tmpAB(3,:) !a-posteriori correction
	tmpE = E
	tmpEx = Ex
	tmpDxx = Dxx
	call DerivativesRHS(tmpAB,tmpE,tmpDxx,tmpEx,tmp)
	do i=1,4
		p2(i,:) = D(i)*(tmpDxx(i,:)+tmpEx(i,:))
	end do
	tmpAB = AB+p2*deltat/2.d0
	tmpAB(1,:) = tmpAB(2,:) + tmpAB(4,:) - tmpAB(3,:) !a-posteriori correction
	call DerivativesRHS(tmpAB,tmpE,tmpDxx,tmpEx,tmp)
	do i=1,4
		p3(i,:) = D(i)*(tmpDxx(i,:)+tmpEx(i,:))
	end do
	tmpAB = AB+p3*deltat
	tmpAB(1,:) = tmpAB(2,:) + tmpAB(4,:) - tmpAB(3,:) !a-posteriori correction
	call DerivativesRHS(tmpAB,tmpE,tmpDxx,tmpEx,tmp)
	do i=1,4
		p4(i,:) = D(i)*(tmpDxx(i,:)+tmpEx(i,:))
	end do
	AB = AB + deltat/6.d0*(p1+2.d0*p2+2.d0*p3+p4)
	AB(1,:) = AB(2,:) + AB(4,:) - AB(3,:) !a-posteriori correction
!	where (AB < 0.d0) AB = 0.d0
return
end subroutine RK
!------------------------------------------------------------------------
subroutine Euler(AB,Dxx,Ex)
double precision, dimension(1:4,1:Ns),intent(inout) :: AB
double precision, dimension(1:4,1:Ns), intent(in) :: Dxx,Ex
integer :: i
	do i=1,4
		AB(i,:) = AB(i,:) + deltat*D(i)*(Dxx(i,:)+Ex(i,:))
	end do
return
end subroutine Euler
!----------------------------------------------------------------------
subroutine erasefiles
	open(unit=1,file=filename0)
	close(1,status='DELETE')
	open(unit=1,file=filename1)
	close(1,status='DELETE')
	open(unit=1,file=filename2)
	close(1,status='DELETE')
	open(unit=1,file=filename3)
	close(1,status='DELETE')
	open(unit=1,file=filename4)
	close(1,status='DELETE')
	open(unit=1,file=filename4b)
	close(1,status='DELETE')
	open(unit=1,file=filename5)
	close(1,status='DELETE')
	open(unit=1,file=filename5b)
	close(1,status='DELETE')
	open(unit=1,file=filename6)
	close(1,status='DELETE')
	open(unit=1,file=filename7)
	close(1,status='DELETE')
	open(unit=1,file=filename8)
	close(1,status='DELETE')
	open(unit=1,file=filename8c)
	close(1,status='DELETE')
	open(unit=1,file=filename8b)
	close(1,status='DELETE')
	open(unit=1,file=filename9)
	close(1,status='DELETE')
	open(unit=1,file=filename10)
	close(1,status='DELETE')
	open(unit=1,file=filename11)
	close(1,status='DELETE')
return
end subroutine erasefiles
!----------------------------------------------------------------------
subroutine savedata(C,R,XI,E,AB,k,time)
double precision, dimension(1:Ns), intent(in) :: C,R,E
double precision, intent(in) :: XI
double precision, dimension(1:4,1:Ns), intent(in) :: AB
integer, intent(in) :: k
double precision, intent(in) :: time
real :: Nlrsp
integer :: i,j,n,front
double precision :: minE
	open(unit=1,file=filename0,position='APPEND')
	write(1,'("Number of save: ",i5,"     Iteration: ",i12,"    Time [min]: ",f12.2)') k,(k-1)*Nsavedata,time*realtime
	close(1)
	n=k
	if (n >= nsave) n=nsave !if more than nsave steps saved, data number nsave overwritten
	!temporary saving
	tmpC(n,:) = real(C(:))
	tmpR(n,:) = real(R(:)/deltat)
	tmpE(n,:) = real(E(:))
	tmpAP(n,:) = real(AB(1,:))
	tmpAM(n,:) = real(AB(2,:))
	tmpBP(n,:) = real(AB(3,:))
	tmpBM(n,:) = real(AB(4,:))
	tmpQ(n,:) = real(AB(1,:)-AB(2,:)+AB(3,:)-AB(4,:))
	tmpAB(n,:) = real(AB(1,:)+AB(2,:)+AB(3,:)+AB(4,:))
	!saving in columns for GnuPlot
	Nlrsp = real(Nl)*real(realspace)
	open(unit=1,file=filename1)
	do i=1,Ns
		write(1,*) (-Nlrsp+real(dble(i)*realspace),tmpC(j,i),j=1,n)
	end do
	close(1)
	open(unit=1,file=filename2)
	do i=1,Ns
		write(1,*) (-Nlrsp+real(dble(i)*realspace),tmpR(j,i),j=1,n)
	end do
	close(1)
	open(unit=1,file=filename3,position='APPEND')
	if (k>1) write(1,*) real(time*realtime),XI
	close(1)
	open(unit=1,file=filename4)
	do i=1,Ns
		write(1,*) (-Nlrsp+real(dble(i)*realspace),tmpE(j,i),j=1,n)
	end do
	close(1)
	open(unit=1,file=filename4b,position='APPEND')
	if (k>1) write(1,*) real(time*realtime),minval(E(:))
	close(1)
	open(unit=1,file=filename5)
	do i=1,Ns
		write(1,*) (-Nlrsp+real(dble(i)*realspace),tmpAP(j,i),j=1,n)
	end do
	close(1)
	!position of the front defined by the minimum of the E field
	minE=10.d0
	do i=1,Nr
		if (E(i) < minE) then
			minE = E(i+Nl)
			front = i
		end if
	end do
	open(unit=1,file=filename5b,position='APPEND')
	if (k>1) write(1,*) real(time*realtime),AB(1,Nl+front)
	close(1)
	open(unit=1,file=filename6)
	do i=1,Ns
		write(1,*) (-Nlrsp+real(dble(i)*realspace),tmpAM(j,i),j=1,n)
	end do
	close(1)
	open(unit=1,file=filename7)
	do i=1,Ns
		write(1,*) (-Nlrsp+real(dble(i)*realspace),tmpBP(j,i),j=1,n)
	end do
	close(1)
	open(unit=1,file=filename8)
	do i=1,Ns
		write(1,*) (-Nlrsp+real(dble(i)*realspace),tmpBM(j,i),j=1,n)
	end do
	close(1)
	open(unit=1,file=filename8c,position='APPEND')
	if (k>1) write(1,*) real(time*realtime),AB(4,Nl+front)
	close(1)
	open(unit=1,file=filename8b)
	do i=1,Ns
		write(1,*) (-Nlrsp+real(dble(i)*realspace),tmpAB(j,i),j=1,n)
	end do
	close(1)
	open(unit=1,file=filename9)
	do i=1,Ns
		write(1,*) (-Nlrsp+real(dble(i)*realspace),tmpQ(j,i),j=1,n)
	end do
	close(1)
	open(unit=1,file=filename10,position='APPEND')
	write(1,*) real(time*realtime),sum(AB(1,:))-sum(AB(2,:))+sum(AB(3,:))-sum(AB(4,:))
	close(1)
	!position of the C front (just for the second graph: position as a function of time)
	front=0
	findfront: do
		if (C(Nl+front) < Cfront) exit findfront
		front=front+1
	end do findfront
	open(unit=1,file=filename11,position='APPEND')
	if (k>1) write(1,*) real(time*realtime),real(dble(front)*realspace)
	close(1)
return
end subroutine savedata
!----------------------------------------------------------------------
end module fonctions

!----------------------------------------------------------------------
!----------------------------------------------------------------------

program liese
use fonctions
implicit none
double precision, dimension(1:4,1:Ns) :: AB !A+,A-,B+,B-,C
double precision, dimension(1:4,1:Ns) :: Dxx !Space derivatives of A+,A-,B+,B-
double precision, dimension(1:4,1:Ns) :: Ex !Derivatives of term with elect. field for A+,A-,B+,B-
double precision, dimension(1:Ns) :: E,R,C !electrical field E, reaction term, C
double precision :: time,XI !time t, current I
integer :: i,j,k
!Initial conditions: square profile
AB(1:2,1:Nl)=a0
AB(3:4,1:Nl)=0.d0
AB(1:2,Nl+1:Ns)=0.d0
AB(3:4,Nl+1:Ns)=b0
C=0.d0
R=0.d0
Ex=0.d0
Dxx=0.d0
call DerivativesRHS(AB,E,Ex,Dxx,XI)
time=0.d0
i=0
k=1
pouet=0
call erasefiles
mainloop: do
	i=i+1
	time=time+deltat
!	call Euler(AB,Dxx,Ex)
	call RK(AB,E,Dxx,Ex)
	do j=1,Ns
		if (infiniterate .eqv. .true.) then
			R(j) = min(AB(2,j),AB(3,j)) !for infinite reaction rate
		else
			if (i==1 .and. k==1) R(j) = rate*deltat*AB(2,j)*AB(3,j)
			if (R(j) > min(AB(2,j),AB(3,j))) R(j) = min(AB(2,j),AB(3,j)) !for finite reaction rate
		end if
		AB(2,j) = AB(2,j) - R(j)
		AB(3,j) = AB(3,j) - R(j)
		if (infiniterate .eqv. .false.) R(j) = rate*deltat*AB(2,j)*AB(3,j) !for finite reaction rate
	end do
	C = C + R
	call DerivativesRHS(AB,E,Dxx,Ex,XI)
	if (i==Nsavedata .or. (i==1 .and. k==1)) then
		call savedata(C,R,XI,E,AB,k,time)
		if (k /= 1) i=0
		k=k+1
	end if
end do mainloop
end program liese
